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We present an extension of the well-known Bogoliubov theory to treat low dimensional degenerate 
Bose gases in the limit of weak interactions and low density fluctuations. We use a density-phase 
representation and show that a precise definition of the phase operator requires a space discretisa- 
tion in cells of size I. We perform a systematic expansion of the Hamiltonian in terms of two small 
parameters, the relative density fluctuations inside a cell and the phase change over a cell. The 
resulting macroscopic observables can be computed in one, two and three dimensions with no ul- 
traviolet or infrared divergence. Furthermore this approach exactly matches Bogoliubov's approach 
when there is a true condensate. We give the resulting expressions for the equation of state of the 
gas, the ground state energy, the first order and second order correlations functions of the field. 
Explicit calculations are done for homogeneous systems. 

PACS numbers: 03.75.Fi,42.50.G 

Recent progresses in the realization of low-dimensional Bose gases in the quantum degenerate regime offer new 
perspectives for the comparison with theoretical treatments. In atomic Bose gases, low-dimensional systems are 
achieved by creating anisotropic trapping potentials. Bose-Einstein condensates of reduced dimensionality, that is 
with the atomic motion frozen in the harmonic oscillator ground state along one or two directions, have been produced 
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Low-dimensional Bose gases have been the subject of early theoretical studies. In the thermodynamical limit for 
spatially homogeneous systems, the Mermin-Wagner-Hohenberg theorem 0, 0| excludes the formation of a Bose- 
Einstein condensate at finite temperature. This is physically due to large phase fluctuations which restrict the 
coherence length of the bosonic field to a finite value. One expects however that strong enough repulsive interactions 
between the particles strongly reduces the density fluctuations of the gas in contrast to the non-interacting case 
1^1^. In this context, Popov introduced a long time ago the concept of quasi- condensate 0- This concept has been 
extended to trapped gases H, 0, . The recent observation of large phase fluctuations for a degenerate Bose gas in 
a highly anisotropic cigar shap ed trap has brought a qualitative experimental confirmation of the theory in a quasi 
one-dimcnsional geometry [llL . 

It turns out that the theory of quasi-condensates has not reached yet the maturity of the theory for condensates. 
In the case of regular 3D Bose-Einstein condensation in the weakly interacting regime, the Bogoliubov theory |l3j |. 
based on a systematic expansion in a small parameter, gives a very precise description of the state of the gas. The 
intuitive idea of the Bogoliubov theory is to use the existence of a single macroscopically occupied mode of the 
field, the mode of the condensate. We recall here the ?7(l)-symmetry preserving version of the theory 0,^3- One 
first splits the bosonic field operator as '0 = 00^0 + ^"05 where uq annihilates a particle in the condensate mode and 
Sip accounts for quantum and thermal fluctuations in the other modes. Then one uses the assumption \dijj\ <^ |ao| to 
solve perturbatively the field equations of motion. This approach is not suitable for a quasi-condensate as there is no 
single macroscopically occupied field mode. Fortunately, in the case of weak density fluctuations, the Bogoliubov idea 
can still be adapted in a quantum phase-density representation of the field operator. One writes the field operator 
^ as exp(i^)p^/^ where 6 and p are position dependent operators giving the phase and the density. One then splits 
the operator giving the density as po + 5p, where po is a c-number and 5p are fluctuations, and one uses the fact that 
/\ ' \5p\ <^ Pq. This idea has already been used in the literature ^] but to our knowledge without a precise definition of 
the phase operator, a well-known delicate point of quantum field theory As a consequence of the non-rigorous 

definition of the phase, divergences appear in the theory 0: one has to introduce an arbitrary energy cut-off, so 
that predictions in ID at zero temperature are made within a logarithmic accuracy only, and in 3D there is no full 
equivalence with the Bogoliubov theory. Another approach based on the current-density operator rather than on 
the phase operator is given by Schwartz |l9j | : an expansion of the Hamiltonian in terms of weak density and current 
fluctuations is performed relating the correlation function of the field to the static structure factor. It is subject to 
the same divergence problem in 2D and 3D in the absence of energy cut-off if one calculates the structure factor in 
the Bogoliubov approximation. 

A possibility to circumvent these difficulties is to rely on the path integral formulation of quantum field theory, 
which involves a functional integral over a classical field, for which the phase is perfectly well defined. This is the 
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approach used by Popov, but with the introduction of an energy cut-off much smaller than the chemical potential 
of the gas, so that the physics at length scales smaller than the healing length is not accurately described. This 
functional integral approach has recently been revisited |^ : it can lead to a cut-off independent formalism for 
quasi-condensates. The predictions of ^2d, ^2^] have the unsatisfactory feature of not reproducing exactly the same 
results as the Bogoliubov theory for a 3D condensate, but this has been corrected in an erratum 22]. 

In this paper, we propose a new improved Bogoliubov approach to treat quasi-condensates in the phase-density 
formalism for a weakly-interacting Bose gas. This approach is based on a lattice model, that is with discrete spatial 
modes, which allows to give a careful definition of the phase operator of the field and to introduce from the start an 
energy cut-off. It uses a systematic expansion in powers of the density fluctuations and of the spatial phase gradient 
and leads to simple expressions for the flrst order and the second order spatial correlation functions of the bosonic 
field that do not depend on the energy cut-off and that exactly reproduces in 3D the predictions of the Bogoliubov 
theory. We also use this formalism to determine the equation of state of the gas to the lowest non- vanishing order in 
the thermal and quantum excitations. 

In section ^ we construct a discretized space model in order to define in a precise way the operators giving the 
phase and the density. We give the physical implications of the space discretization, restricting this approach to 
highly degenerate and weakly interacting Bose systems. In section ^ we derive a quadratic approximation to the 
Hamiltonian, that is we derive approximate linear equations of motion for the density fluctuations and the phase 
operators. We recover to the lowest order the Gross-Pitaevskii equation for the quasi- condensate density and we 
recover the Bogoliubov spectrum for the excitations. We also push the expansion to the next order, by producing 
a cubic correction to the quadratic Hamiltonian, including the interaction between the quasi- condensate and the 
excitations. We show that inclusion of this correction is necessary to get a consistent theory and to establish the full 
equivalence between our approach and the number conserving Bogoliubov theory. In section IIIII we present a few 
applications of our formalism: we give general formulas for the equation of state and the ground state of the gas, and 
for the flrst order and second order correlation functions gi and 52 of the field operator. In section Hvl we apply our 
formal results to the homogeneous Bose gas in various dimensions of space. This allows to derive simply the validity 
condition of the method and to compare with existing results in the literature. 

I. CONSTRUCTION OF A DISCRETE PHASE-DENSITY REPRESENTATION 

A. Why discretize the real space ? 

In previous studies of quasi-condensates the basic tools of the theory are an operator 

p{v) = i^Hv)ij{v) (1) 

giving the density in r and an operator ^(r) giving the phase of V'(r), the field operator in r, the position r being a 
continuous variable p^ . A small parameter of the theory characterizing the regime of quasi-condensates is then that 
the density fiuctuations, that is the fiuctuations of /5(r), are small in relative values: 

Var(p(r)) ^ {p{vf) - {p{v)f « {p{v))\ (2) 

However one finds that the expectation value of /5(r)^ is infinite in every point with a non-vanishing mean density 
p(r) = (p(r)): 

{p{vf) = 5{Q)p{v) + (^t(r)^t(i.)^(r)v}(r)) (3) 

where (5(0), the value of the Dirac distribution at the origin, is infinite, and the second term in the right hand 
side, giving the probability density of finding two atoms in the same point of space, is finite in any realistic model. 
Mathematically this divergence is due to the use of the bosonic commutation relations of the field operators ip{v) and 
■!/;^(r) in the same point of space to put the atomic field product in normal order. 

In order to have small, and therefore finite, density fluctuations, one is forced to discretize the space, that is to 
collect the particles in little boxes at the nodes of a spatial grid. Each little box has equal length I along each dimension 
of space and is parameterized by the position r of its center. The held operator '0(r) has the effect of removing a 
particle in the box at position r and it now satisfles the bosonic commutation relations: 



[^(r),^t(r')] = %: 



(4) 
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where dry is the discrete Kronecker delta and D is the dimension of space. The variance of the operator giving the 
density is now finite: 



Var(p(r)) = {^^^ {r)i^\r)i^{r)i^{r)) - p^r) 

ID ■ 



(5) 



In the vahdity domain of the theoretical approach of this paper, this variance will be much smaller than p^{r) because 
both the term on the first line of the right hand side and the term on the second line are small: 

|(^t(r)v;t(r)V'(r)^(r)) - p'{r)\ « p\r) (6) 

p(r)/^ » 1. (7) 



B. The phase operator 

In usual continuous space theories an hermitian field phase operator 9{r) is introduced subject to the following 
commutation relation with the operator giving the density: 

[p(r),^(r')]-*<5(r-r'). (8) 
In our discrete model the desired commutation relation is modified into 

[pir),§ir')]=^^-jf. (9) 

First we recall briefly that there exists actually no hermitian operator 0{r) satisfying strictly the above commutation 
relation. From the identity ^ one can indeed show that the operator 

r(a) = e-^"^^'"), (10) 
where a is any real number, is a translation operator for the density : 

T(a)tp(r)r(a)=/5(r) + ^. (11) 

This identity contradicts two fundamental properties of p{r), the positiveness and the discreteness of its spectrum 

m 

We now proceed with the construction of a phase operator 9{r) approximately satisfying the commutation relation 
0. The key ingredients allowing such an approximate construction are (i) to be in the limit of a large occupation 
number of the considered box of the lattice, and (ii) to construct the operator e'^ first, which, according to (|ll|l taken 
with a — —1, simply reduces the number of particles in the considered box by one. 

In each spatial box we introduce the basis of Fock states |n, r) with exactly n particles in the box. In this basis the 
field operators have the following matrix elements: 

V>(r')|n,r) ^/n\n - l,r), 

4^\r')\n,r) = -iil V^T+T |n + 1, r) (12) 

as a consequence of the commutation relation Q). The atomic density p defined by p{r) — ip'' {r)'ip{r) is diagonal in 
the Fock state basis: 



Tl 

/5(r')|ri,r) = 5r,r' To |f^,r). (13) 



We then introduce the operator A defined by: 



V'(r) EE Air)^^. (14) 
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In the Fock space A{r) reduces by one the number of particles n in the box r: 

A{r')\n,r) = [l - Sr,fl) Sr^r'\n - l,r). (15) 

Note that its action on the vacuum state of the box gives zero. For each box r, the definition of A leads to the exact 
relations: 

AA^=I, iU = /-|0)(0| and = |0)(0|, (16) 

where / is the identity operator and |0) is the zero-particle state or vacuum state in the box of center r. We find 
that the operator A is almost unitary, i. e. it is effectively unitary for a physical state of the system with a negligible 
probability of having an empty box. In what follows, we assume that this condition is satisfied, so that the projector 
|0)(0| can be neglected: 

occupation probability of \n = 0, r) ^ 1. (17) 
In this case, we write the approximately unitary operator A as: 

i(r)~e^^W with ^t(i.)~^(r) (18) 
which amounts to writing the field operator as: 

^(r)~e'«»y;5(0. (19) 

This should be understood as a formal writing, allowing for example to recover the matrix elements of A and therefore 
of the field operator from the commutation relation @ . We summarize below all the commutation relations of our 
phase-density representation: 

[p(r), eir')] ^ ^ [p(r), /5(r')] = [^(r), 0(r')] (20) 

We come back to the constraint (|17|l at the basis of the construction of exp{i9). A sufficient condition to have a 
low probability for a zero particle occupation in a box is obtained for a large mean number of particles in the box and 
with small relative particle number fiuctuations. This is the regime that we wish to consider in this paper. We are 
therefore back to the discussion of the previous subsection and to the conditions (|6I7|) for weak density fiuctuations. In 
particular, the construction of the operator exp(i6') becomes problematic in the limit I — > 0, that is in the continuous 
model. 



C. How to choose the grid spacing I 

Working on a grid can be also seen as performing a coarse-grain average over all physical quantities on a scale I. 
This averaging suppresses the short wavelength modes (shorter than I) and thus introduces an energy cutoff: 

Ecut^-^. (21) 
ml 

This cutoff is of no physical consequence if all characteristic energies (/i, ksT) are smaller, i.e., that I is smaller than 
the corresponding characteristic lengths. This leads, for instance, to the following restrictions for I: 

and 1<X (22) 

where 

(23) 



is the healing length, and 



L. T (24) 
mkB 1 
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is the thermal de Broghe wavelength. These two restrictions, combined with Eq. impose: 

pX" > 1, (25) 
» 1- (26) 

These are conditions of validity for our discrete model. 

The first one, Eq. I|25|) . is the quantum degeneracy regime occurring at sufficiently low temperatures. The second 
restriction, Eq. (|26|l . corresponds to the regime of weakly interacting systems. Its dependence on the density varies 
according to the dimension of space. In ID and 3D, the mean field prediction for the chemical potential is /i ~ gp where 
(7 is a constant characterizing the interaction potential, the so-called coupling constant. In ID, Ea. (|26|l is the high 
density limit where a mean field theory is valid; we recall that the small density limit p^ = fi^J p/mg <C 1 corresponds 
to the strongly interacting (or strongly correlated) Tonks gas regime. In 3D, the effective coupling constant g is related 
to the s-wave scattering length a of the interaction potential, g = Anh'^a/m, so that cx pa^ ^ 1: one recovers 
the usual small gaseous parameter \fpcP' . In 2D, the chemical potential scales as ?i^p/(mln(l//9a^)) where a is the 
scattering length of the 2D interaction potential so that the condition p^ ^ 1 results in a low density condition, 
Hl/pa^) > 1. 



II. PERTURBATIVE TREATMENT OF A MODEL HAMILTONIAN 



A. Model Hamiltonian 



In our lattice model, we represent the binary interaction potential among the particles by a discrete delta potential: 



V{ri - ra) = ^(^n.r^ 



(27) 



where go is the bare coupling constant. Note that go in general differs from the effective coupling constant g, and we 
shall come back to this point in section UlI Al With this model potential, the grand canonical Hamiltonian is: 



^ = E [ - ^i'Hr)A^j{r) + ([/(r) - p)iJ^\rUj{r) + |^t(r)^t(r)^(r)^(r) 



(28) 



where U{r) is an external trapping potential and where the Laplacian is a symmetric operator coupling the different 
neighboring boxes: 



A/(r)=5] 



/(r + Zej) + /(r-?ej)-2/(r) 



(29) 



The Gj arc the unitary vectors and j the different orthogonal space directions (for example, j — x,y,z in 3D). As 
usual we take periodic boundary conditions inside a rectangular box with lengths integer multiple of I. 

We now rewrite the Hamiltonian in the density-phase representation, that is in term of the operators p and 9 giving 
the density and the phase as defined in the previous section. The contributions of the trapping potential and of the 
interaction potential to the Hamiltonian are local in real space and therefore include the operator p only: 



C/(r)-p + |U(r)-^ 



(30) 



where we have used the bosonic commutation relation Q to exchange one of the with in the interaction term. 
The kinetic energy term involves explicitly the phase operator: 



kin 



2ml'- 



(31) 



where we have introduced the notations 9±j = 6{r± Zej) and p±j ~ /o(r ± ^ej). A remarkable property of this writing, 
to be used below, is that it involves only the difference of two phase operators between two neighboring points of the 
lattice. 
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B. Hamiltonian quadratisation and cubisation 



We now expand the Hamiltonian to third order in powers of two smah parameters. As aheady discussed in 
subsection II Al the regime of quasi-condensates that we are interested in corresponds to smah relative fluctuations 
5p of the density. In the zeroth order approach totally neglecting the density fluctuations, the density is set to a 
deterministic value poi as we shall see. The second order expansion allows us to describe the density fluctuations: 

p(r)=po(r)+(5p(r). (32) 

The third order expansion allows us to calculate the mean value of 5p{v). 

The first small parameter of the systematic expansion used in this paper is therefore given by 

1. (33) 

where \5p\ is the typical value of the operator 5p in the physical state of the system. Mathematically this allows to 
expand Vp as 

"1/2 1/2 ^ 1 ^ 5p^ I Sp^ 

The second small parameter of the expansion is given by 

£2 = |;v^| < 1. (35) 

Here V represents the gradient on the lattice: 

V;(,)^^ /(- + ^-^)-/(--^'-^) e, (36) 
j 

where / is an arbitrary function. Physically the existence of the small parameter €2 is reasonable: it is at the basis 
of our hypothesis that the continuous quantum field problem can be well approximated by a discrete lattice model, 
provided that I is small enough, see subsection II CI Mathematically this second small parameter allows to expand in 
(|3I|I the exponentials of the phase differences: 

e'^'^^-'^ - 1 + ^ih, -9)- l{9+j ~6f... (37) 

From the fact that the discretisation length / is on the order of the smallest of the two macroscopic length scales 
^ and A, see (|22f) . it will be checked later that the parameters ei and e2, though of an apparently different physical 
origin, can be chosen to be of the same order of magnitude, 

,38, 

and can therefore be treated mathematically as infinitesimals of the same order. The mathematical details of the 
expansion 



H = Ho + Hi+ H2 + H:i + 



(39) 
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are given in the appendix^ we present here only the results: 



Hn = 



Hi = 



E [ - ^^^^ + f ^0 + (Uir) - m)Po^ 

r 



^y/po + U{r) -11 + goPo Sp 



H2 = E2[po]+Y,l''[- 



Sp 



2m 2^ V2VA) 



jj^A^ + -5p 

ompQ ^ 



— E V^o(r)po(r + /ej; 



P 



Aral'- 



1/2 



1/2 



Po Po^+j J 



n X^Op ( -1/2. 1/2 1/2. -1/2 

+ ^E-(Po APo -Po APo 



16m 



E^" 



5p^ 

5/2 

LPo 



avp^ 



pT 



A 



^P 



(40) 



The quantity E2 in i?2 is a c-number functional of the density poi given in the appendix^ which has therefore no 
contribution to the dynamics of the quantum field. 



C. Iterative solution for the quadratic Hamiltonian 



We now solve perturbatively, order by order, the Hamiltonian problems defined by iJp, ^^o + Hi and Hq + Hi + H2. 
To zeroth order in £1^2, the Hamiltonian is a c-number. As the chemical potential is fixed in our approach, Hq is 
minimized for a density profile po(r) such that ^fp^ solves the discrete version of the Gross-Pitaevskii equation: 



'7^A + U{v) - p + gapa 
Am 



This density profile constitutes the zeroth order approximation to the density p. 
that we call iVg: 



(41) 



ft contains a number of particles 



N,^Y.^^p,{v). 



(42) 



Note that A'o coincides with the mean total number of particles N only to lowest order in the theory. Equation H41|l 
defines po and therefore iVo as function of the chemical potential p. It will however turn out to be more convenient 
to parameterize the theory in terms of A'o rather than in terms of p. We will therefore consider p and po a-s functions 
of A^o: 



Po(i-) 



Aio(^o) 
Po(r;A^o)- 



(43) 
(44) 



po is therefore the Gross-Pitaevskii prediction for the chemical potential of a gas of Nq particles. 

For the choice of density profile H41(l . the first order correction Hi to the Hamiltonian vanishes. We therefore have 
now to solve the Hamiltonian problem defined by H2, in order to determine the lowest order approximation to the 
density fluctuation 5p and the phase 9. It is instructive to write the corresponding Heisenberg equations of motion, 
which are linear (and therefore trivially solvable) since H2 is quadratic. As 9 and 5p are two canonically conjugate 
variables, the equations of motion are: 



ndt9 
hdtSp{r) 



1 dH2 
ID dSp{r) 
1 dH2 
d9{T) 



2mJ^ 



A 



5p 



A ie^o 



-Sp 



AJp^ 



2po 
9A^ 



- 9o5p 



(45) 
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An important difference of these equations with the so-called quantum hydrodynamics equations for 5p and 9 is that 
our formalism keeps the so-called quantum pressure term for dt9 whereas it is usually neglected in the literature |23j |. 
This allows our treatment to have a cut-off energy larger than /i, whereas the usual treatment is restricted to energy 
modes much below /i. 

Furthermore one can simplify these equations using the Gross-Pitaevskii equation l|41|l to eliminate Ay^pg: 



fidtO 
hdtSp{r) 



1 



- — A + U + 3gopo - M 
2m 



Sp 



2m 



A + U + gopo - ^l 



(46) 
(47) 



This gives the idea of a very simple canonical transformation which, remarkably, maps our equations for a quasi- 
condensate H46I47|I into the equations for the Bogoliubov modes of a condensate: the field 



B 



Sp 
2Jp^ 



+ Wpo 



(48) 



has bosonic commutation relations 



[B{r),B^r')] 



ID 



and it obeys the standard Bogoliubov equations 



ifidf 



B 

B^ 



B 
fit 



I^A + U-^l + 2.goPo 



2m 



.90 Po 

A + [/ - /i + 2goPa 



B 



(49) 



(50) 



This mapping can be readily extended to the Hamiltonian H2, which is expected to be canonically equivalent to the 
Bogoliubov Hamiltonian: 



Ho 



idj2b^ 



2m 



A + U + gopo ~ p] B + gopo 



B^B + ^ {B^+B^'^ 



(51) 



We have checked that the identity H51I) indeed holds by replacing B by its expression (|48l) in terms oi Sp and 6, and 
by using the value of the commutators (|2()|l and the fact that solves the Gross-Pitaevskii equation. Remarkably 
the energy functional E2[po] is exactly compensated by the contribution of the commutators. 

This mapping therefore allows to reuse the standard diagonalisation of the Bogoliubov Hamiltonian. We recall here 
briefly the procedure described in n^l27j. One introduces the normal eigenmodes {us, Vs) of the Bogoliubov operator 
£gp with an energy e^, normalisable as 



^/^[|«.(r)p-K(r)p] = l. 



(52) 



Then (w*,m*) is an eigenmode of £gp with the energy — e^. To form a complete family of vectors one has to further 
introduce the zero energy mode of £gp, given by ((j)Q, ^(f>o), and the anomalous mode (0a, (j)a) with 



00 = Vpo/No 



and 



^dNo\/PO 



(53) 



The corresponding normalization of the anomalous mode is such that the scalar product of 0o and 0q is 1/2. With 
these definitions, one introduces the components of {B, W) on the zero energy mode, on the anomalous mode and on 
the regular {us,Vs) modes: 



P 

7K 



-bi 



(54) 



Q is a collective coordinate representing the quantum phase of the field and P is its conjugate momentum 

[P,Q] = -i. (55) 
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Physically P corresponds to fluctuations in the total number of particles, as expected, and as shown in more details 
later, see 1)6 7|l . The operators 6s are bosonic annihilation operators with the usual commutation relations [6^,6^'] = 
(5s. s'. They commute with P and Q. The inverse formulas giving bg, Q and P in terms of B can be found for example 
in |l5| . Equation l|54|l results in the following modal expansion for the density fluctuations and the phase operators: 

^(r) = J2^,{r)bs + eUr)bl-Q 
5p{v) - ^<5ps(r)Ss + <5p:(r)6t+pa^^po ^^^^ 



where 

Os{v) 



Usjr) - Vsjr) 
2VPo(r) 



Sp,{r) = v/po(r) (usir) + v,{r)). (57) 
By construction, this modal expansion, when inserted into the quadratic Hamiltonian H2, results in 

H = Y^ ^sb% + ipVo + E2[po], (58) 

s 

where /ig = dfio/dNo- This is the sum of uncoupled harmonic oscillators, plus a massive free degree of freedom 
corresponding to the unbound phase variable Q. The effective mass of the phase variable is given by I/pq. The 
energy functional i?2[po] will be calculated in llll Bl where it will be shown that it leads to exactly the same ground 
state energy as the number conserving Bogoliubov theory. This shows that the Bogoliubov theory can be used to 
calculate the ground state energy even for e.g. ID quasi-condensates, a fact commonly used in the literature [2^ |29| 
but which looks rather heuristic in the absence of justiflcation. 



D. Effect of cubic Hamiltonian corrections on the density 



The physics contained in the cubic term H3 of the Hamiltonian is very rich. It includes interaction effects between 
the Bogoliubov modes of the previous section, allowing a generalization to quasi-condensates of the theory of energy 
shifts and Beliaev-Landau damping usually put forward for Bose-Einstein condensates [s^, |^ 112 • 

We are more modest here. Our motivation to include the cubic corrections is that the quadratic Hamiltonian H2 
brings actually no correction to the zeroth order approximation po to the mean density, since the mean value of Sp 
vanishes at the level of the second order theory. This is highly non satisfactory as it brings some inconsistency in the 
calculation of an observable like gi, the first order correlation function of the fleld: to get a non-trivial prediction for 
gi one has to include terms quadratic in the phase operator, which are second order in €2, which forces to also include 
second order corrections to the mean density, as will become very explicit in section UTTI 

We therefore calculate the first order correction to the equations of motion of 5p and 9 due to the cubic Hamiltonian 
term H3, and we take the average over the quantum state corresponding to the density operator at thermal equilibrium 
for the Hamiltonian H2. This gives source terms to add to the equations for the mean density and phase derived from 
H2. We leave the details of the calculations for the appendix 151 and we give directly the result: 



fidt{5p)z = 



1/2 

Po 



2m 



A + U + gopo - p 



1/2 



U + SgoPo - p 



{^pl^'iO), 
{Sp)3-{B^B)2^ 



1/2 
Po 



+ gopri^B^B + B^ + St2)2 - 2{P^)2p'„dNo 



(59) 



(60) 



where the thermal average (. . .)2 is taken with the unperturbed Hamiltonian H2 and (. . .)3 is taken with the perturbed 
Hamiltonian H2 +H3 to first order in H3. The expectation value of the 'kinetic energy' of the unbound phase variable 
in H58() is equal to kBT/2 according to the equipartition theorem so that 



Mo 



(61) 
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At equilibrium the expectation value of dtSp and dt9 vanish. This fact is obvious for dt5p\ it is less obvious for 
dtO because of the presence of the unbound variable Q, we therefore produce a proof of that in the appendix |0 We 
therefore have to solve (|59I60|I with the left-hand side set to zero. The first equation H59|) imposes that the mean value 
of 9 is position independent, a trivial result. In the second equation, the operator acting on {5p)^ is strictly positive 
so that it is invertible and H60(l determines the correction to the mean density in a unique way. 

We now go through a sequence of transformations allowing us to get a physical understanding of the value of {5p)^. 
The first step is to pull out the contribution of the 'anomalous' terms P, Q in the modal expansion (|54|l : 



(62) 



We calculate the expectation values of H60|l involving the operator B, using the fact that all the crossed terms between 
the anomalous part and the operators 6s have a vanishing expectation value: 



{B% + (i3t2) 



Nn 



(63) 
(64) 



The term (t>Q(t>a in H63|) comes from the non-commutation of P and Q, see (|55|) . The contributions of {0^)2 in H63|l 
and H64() . when inserted into H60|l . are shown to compensate exactly when one uses the fact that 00 solves the Gross- 
Pitaevskii equation. This was expected from the C/(l) symmetry of the Hamiltonian: only differences of the phase 
operator in two points appear in the Hamiltonian, so that H does not depend on Q and the mean density does not 
depend on {Q^)2- 

We therefore get an equation for {Sp)^ involving the expectation value of as a source term, and which looks 
rather involved: 



= 



h 

-—A + U + 3gopo - p 
Zm 



'{Sph - 4'lNo\P^)2 - {{BlBn)2 - Ma)^ 



1/2 

Po 



gapl'^{A{BiB,, - + K + bI^ + {P% 6gop'J^dN„^~ 2p[, d^oVP^ 



(65) 



Fortunately the underlying physics is very simple and allows to predict the effect of this source term on the mean 
density. One first identifies the physical meaning of P in 1)56(1 . Using the well-known fact that the eigenmodes of £gp 
are orthogonal for the modified scalar product of signature (1,-1), one has 15]: 



t>o\Us 



t>0\Vs 



^;^0o(r) K(r)+z;,(r)] =0 



(66) 



so that the sum of 6ps over all spatial nodes vanishes. As a consequence the operator N giving the total number of 
particles in the gas is simply 



N: 



^Z^p(r) = iVo + P 



(67) 



where we have used the identity 



5:;^aA.„po(r) = ^$:/^Po(r 



dNo 
dNo 



(68) 



The source terms involving P therefore correspond to fluctuations in the total number of particles in the gas, due 
to the fact that we consider the grand canonical ensemble. The effect of these grand canonical fluctuations can be 
considered for the case of a pure quasi-condensate at the order of the present calculation so it is easy to calculate it 
directly. In the grand canonical ensemble the probability that the quasi-condensate has n particles is 



n„ cx exp [-P {EQ{n) - pn)] 
where Eo{n) is the Gross-Pitaevskii energy for the density profile po{r;n): 



Po{r■,n)A^y po{r;Ti) + U{r)po{r; n) + ^pl{r-n) 



(69) 



(70) 
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The corresponding mean grand canonical density is 

Pgc(i") = / dnIlnpo{r]n) 



(71) 



where we treat n as a continuous variable. The zeroth order approximation n = Nq for the number of particles in the 
quasi-condensate is such that Eo{n) — fin has a minimum: 



— [Eo{n) — fin] = fJia{n) — /i = for n = Nq 



(72) 



as shown in l|43(l . The corresponding density profile is pQ{r;No). The next order correction to that is obtained by 
expanding the n-dependent density profile to second order in n — Nq and by averaging over n: 



SpGcir) = {{n-No))dN„Po{r;No) + -{{n ~ Nof)d%^^poiv; No). 



(73) 



The second moment of n — A^o is calculated to lowest non- vanishing order by a Gaussian approximation to 



1 E 1 
Eo{n) -pn~ Eo{No) - fJ.No + 3 "^(^ " ^0)^ = const + -fi'oin - Nof 



(74) 



This leads to 



{{n~Nof 



' Gauss 



Mo 



/p2 



(75) 



More care has to be taken in the calculation of the mean of n — No'- the Gaussian approximation to n„ gives a 
vanishing contribution, so that the cubic distortion to it has to be included: 



Eo{n) -fin ^ Eo{No) - fiNo + ^fi'o{n - Nof + i/io(" " ^0)^ 



n„ cx exp 



lpf,'^(n-Nof 



l-i/3Mo(«-^o)' 



(76) 
(77) 



We then get a non- vanishing mean value for n — No 

1 



({n - A^o))distor = "77/3Aio(('^ ^ ^0) ) Gauss 

6 



(78) 



We have therefore predicted in a very simple way the correction to the mean density due to grand canonical fluctua- 
tions: 



1 



d%^,Po{r;No)-^dNoPo{r;No) 
Po 



(79) 



How does this compare to the general formalism (|65|) ? We need to obtain a partial differential equation for SpQQ. 
We just take the second order derivative of the Gross-Pitaevskii equation (|41|l with respect to iVo and we replace po 
by a/po^ in the resulting equation and in H79|) . This leads to the remarkable identity 



2m 



A + U + 'igoPo - p 



5pGC-N^'cf>l{P^)2 
1/2 
Po 



P^)2 Qgopl^^dNoVPO -'^Po dNo^/PO- 



(80) 



The right-hand side of this identity coincides with the source term of (|65|l involving P\ We have therefore successfully 
identifled 5pQc as a piece of (^p)3 and we are left with the simpler equation 



- — A + U + igoPo - P 
Zra 



1/2 

Po 



9qPo 



1/2 



mix - 4>0<t>a) + Bl + i3t2)2. 



(81) 
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We are not totally satisfied yet since the operator i?„ does not obey bosonic commutation relations when the system 
is not spatially homogeneous, in particular the field _B„ does not commute with itself when taken in two different 
points: 



[i3„(r),B„(r')] - 0a(r)0o(r')-0a(r')0o(r). 
[i3„(r),Bt(r')] = ^5.y-Mv)4>a{v')-M^')M^)- 



(82) 
(83) 



To circumvent this difhculty we split the field i?„ in its component along the quasi-condensate mode 0o and its 
orthogonal component: 

^„(r) =a(/.o(r) + A(r). (84) 
The Bogoliubov functions Us{y) and Vs{v) can be chosen here to be real. The operator a can be then written as 



(85) 



where we have used the property (|66|l . This clearly shows that the operator a is anti-Hermitian: 



= —a. 



The field A has the following expansion on the 6^: 

A(r) = u,^{r% + v,^{r)bl 



(86) 



(87) 



where the index _L indicates projection orthogonally to (j)Q. This field has now the desired bosonic commutation 
relations 



[A(r),A(r')] - 

[A(r),At(r')] = ^S^y ~ Mr)Mr'). 



(88) 
(89) 



Note that a does not commute with A: 



[d,A(r)] = -0o(r)-0a(r). 



(90) 



We insert the splitting of _B„ in (|81|l . The terms quadratic in a cancel exactly, in the same way the terms in 
canceled. The terms linear in a can all be expressed in terms of the expectation value of an anticommutator, 
({d,A})2 using the commutation relation (|90|) and the fact that (dA)2 is a real quantity. Furthermore, using the 
techniques of Appendix E of p^ . as shown in the appendix^ one obtains a simple partial differential equation for 
the anticommutator: 



- — A + U + goPo - M 
2 m 



({d,A(r)})2 = -5];^5oPo(r')0o(r')({A(r')+At(r'),A(r)})2. 



(91) 



Remarkably this allows to eliminate completely the operator a in (|81|l ! We finally get an equation for the correction 
to the mean density involving the operator A only: 



= 



-—A + U + Sgopo - fi 
2m 



{6p)3 - SpGc " (AU): 



Sir) 



(92) 



where we have introduced the source term 

Sir) EE .goA^o0o(r)(4At(r)A(r) + Pir) + A^^{r) - ^^(r))^ - ^ .goPo(r')0o(r')({A(r') + At(r'), A(r)})2. (93) 



It will reveal convenient to introduce the function x(r) defined in a unique way by 



= 



- — A + U + 3gopo - p. 
2m 



Xir) + ^Sir). 



(94) 
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We then obtain the following final expression for the correction to the mean density due to the cubic Hamiltonian 
terms H^: 

{Sp)3{r) = SpGcir) + 2^o(r)x(r) + (A+(r)A(r))2. (95) 

In the particular case where the gas is Bose condensed, our general theory for quasi-condensates also applies, of 
course. One then expects that the result (|95() has already been obtained for the condensate and can be given a 
clear physical interpretation. This expectation is totally justified: as shown in the appendix El the component of 
x(r)/-/Vo orthogonal to (t>a is the correction given in to the Gross-Pitaevskii condensate wavefunction due to 
the interaction with the non-condensed particles; the component of x along (j)Q describes the condensate depletion, 
and (A+A)2 is the mean density of non-condensed particles. 



III. APPLICATIONS OF THE FORMALISM: GENERAL FORMULAS 



A. Equation of state 

What is referred to as the equation of state of the gas is the expression of the chemical potential as function of 
the mean total number of particles N and the temperature T. It is useful in particular to predict properties of an 
inhomogeneous gas within the local density approximation. 

We therefore have now to calculate ^ for the quasi-condensate. This is equivalent to a calculation of Nq as and 
A^o are by definition related through 143(1 . To lowest order of the theory one assumes a pure quasi-condensate with a 
density profile p(r) ~ Polr), where ^/po solves the Gross-Pitaevskii equation H41|) . One therefore gets N — Nq so that 

^^ = MN)■ 

The first non- vanishing correction to the density profile is given by 1(9 5|l . By integrating ((95(1 over space we get the 
corresponding correction for the mean total number of particles: 

N = Nq + 6N (96) 

r r 

The contribution to 6N due to our use of the grand canonical ensemble can be calculated exactly from a spatial 
integration of ifT?^. using the same technique as in (|()H|l : 

SNgc = -ksT^. (98) 

The contribution of the term involving x can also be made explicit by multiplication of l(94() by (/'a(r) defined in 1(53(1 
and by spatial integration. The function 0^ (r) is indeed known |27j to solve the partial differential equation 



- — A + U + 3goPo - M 
Zm 



= A^oMo'^o(r) (99) 
which can be checked easily, just by taking the derivative of ((41(1 with respect to A^o- This leads to 



J2 2Mr)xir) = -T7^ E l''Mr)Sir) (100) 

where the source term S is known explicitly, see 1(93(1 . We just have now to replace A^o by A^ — 6N in ((43(1 and expand 
to first order in 6N: 

p = fioiN-SN)^^ioiN)-6Nfi'oiNo). (101) 
We obtain the following expression for /x: 

^ MN) + kuT^ - ^^'oiNo) J2 Z^(Anr)A(r))2 + ^ E ''''^a(r)^(r). (102) 
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Equivalently we can replace the source term by its explicit expression to get 

^ c MN) + ksT^^ - ^i',{N,) + E ^''(At(r)A(r))2^ + ^ g, (9^„po(r)) (2(AtA)2 + Re(A2)2) 

-5]l%</.g(r)({A(r)+At(r),7})2 (103) 

r 

where we have introduced the operator 

7 = ^/^0,(r)A(r) (104) 

r 

and we have used the identity 

Mo = E ^''3o0^(r)ajVoPo(r; N^) (105) 



obtained by performing the scalar product of both sides of H99f) with (/)o. Application to spatially homogeneous systems 
will be given in section Hvl in this case both the operator 7 and /Iq vanish. 



B. Ground state energy 

We now show that the ground state energy of a quasi-condensate can be calculated with exactly the same Bogoliubov 
formula as for the ground state energy of a condensate. 

We have to determine the ground state energy of H2 ■ We write it as the expectation value of 1)51(1 at zero temperature, 
that is here in the vacuum of the 6, and of P: 



grounc 



A{H2) = l^'Y.^B^ {"^^ + ^ + ^"-"o " ^) + 5oPo( B^B + \{b^ + B^^ 



(106) 



We reproduce the transformation of subsection IIIDI We split B in an anomalous part involving P, Q, plus the 
contribution of the antihermitian operator a and of A, the orthogonal component of the normal part. In the first 
expectation value of the right-hand side of (|106|l the operators Q and a disappear as they come with the factor (/'o(r) 
in B and (pQ solves the Gross-Pitaevskii equation (|41|l . The expectation value of P^ in the ground state of H2 also 
vanishes, so that 



(St ( - — A + U + gopo - M j B)2 = (A^ f - — A + U + goPo - M ) A)^ 



(107) 



The same transformation is applied to the last expectation value in ((l()6|l . Remarkably the terms involving a exactly 
cancel when one uses the relations and (|D5p . This leads to 



AU + - ( A2 + At2 



The expectation values involving A are readily calculated from the modal expansion l|87|l : 



ground(iJ2) = - J E '^.9oPo0o + E< 



-9qPq\Us±) 



(108) 



(109) 



As (us, Vs) is an eigenvector of £gp, (usJ_, Vs±) is an eigenvector of the operator £ defined in and this expression 
can be further simplified to 



ground(i72) = "i^^l'^ 9oPo4'l -^es{vs±\vs±). 



(110) 



The last step is to include the contribution of Hq and to remove the —pN term from the grand canonical Hamil- 
tonian. The ground state energy of the canonical Hamiltonian for N particles is therefore 



£^grou„d(iV) :^pN + EoiNo) - pNo + ground(i/2) 



(111) 
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where Eq is the Gross-Pitaevskii energy H70|l . As we did in subsection IIII Al we replace A'o by — 5N, where SN is 
calculated from H^, and we expand Eq{N — 5N) to first order in 6N: 



^lN + So(A^o) - M^^o - Eo{N) - 6Nipo{N) - ^) ~ EoiN). 



(112) 



We recall that by definition fi = p,o{No). The first term in the right hand side of pi()|l amounts to performing a small 
change in the Gross-Pitaevskii energy functional, expressing the fact that a given particle interacts in the gas with 
A^ — 1 particles so that the mean field term should be proportional to A^ — 1 rather than to A^ . The final expression 
for the ground state energy is: 



-|-<^o(r; N)AMr; N) + U{v)4>l{v; N) + ]^g^{N - l)<j,l{v- N) 



r 

This exactly coincides with the Bogoliubov result, see e.g. equation (71) of [l 



- ^es(wsj_bs±). (113) 



C. Second order correlation function 



The second order correlation function of the atomic field is defined as: 

g2(r) EE (^t(r)^t (0)^(0)^(1.)) (114) 

where we have taken for simplicity one of the two points as the origin of the coordinates. To calculate 52 with the 
formalism of this paper we have to express g2 in terms of the operator p giving the density. This is achieved thanks 
to the commutation relation Q of the bosonic field i/;: 

g,ir)^{pir)m)-jn-{m)- (115) 
We then insert the splitting (|32|) of p in terms of the quasi-condensate density po and the fiuctuations Sp: 

52(r) = Po(r)po(0) + Po(0)(<5/5(r)) + po(r)(5p(0)) + (<5p(r)^p(0)) ~ ^ [po(0) + (<5p(0)>] . (116) 

This expression of 32 is still exact. We now perform approximations consistent with an expansion of g2 up to second 
order in the small parameters £1^2- The expectation value of the term quadratic in Sp is calculated within the thermal 
equilibrium for the quadratic Hamiltonian i?2- The expectation value of 5p is evaluated in subsection lll DI bv inclusion 
of the cubic perturbation H3. The contribution of 5p in the last term of (|116|) is negligible as it is ef times smaller 
than the leading term in g2. We therefore obtain the explicit expression 

32(r) ~ po{r)po{0)+PoiO){Sp{r))3+po{v){Sp{0))3 + (<^/5(r)<5p(0))2 - jj^Po{0). (117) 

This writing is however not the optimal one as the last term in l/l^ gives the wrong impression that 52(0) strongly 
depends on the discretisation length I in the continuous limit Z — s- 0. In fact this strong dependence exactly compensates 
a term in l/l^ in the density fluctuations {Sp^{0)) coming from the fact that Sp^{0) is a product of field operators 
not in normal order. To reveal this fact we express Sp in terms of the operator A of (j87|l 

Sp{r) = (A(r) + At(r)) + PdNoPo{r; No) (118) 

and we put the resulting expression in normal order with respect to the field A thanks to the bosonic commutation 
relation l|S^ : 

Sp{r)Sm = ■.Sp{r)Sp{0) : +j^po{0) - A^o0^(r)0^(O) (119) 

where : : is the standard notation to represent normal order. The spurious term in l/l^ is then exactly canceled: 

g^ir) = A^o(A^o - l)0o(r)0o(O) 

+ po{0){Sp{r))3+po{r){Sm)3 + {-Sp{r)5m-)2- (120) 
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This expression allows to prove the equivalence with the prediction for 52 in the Bogoliubov theory. We do not present 
the calculations here, as they are a straightforward application of appendix ^ Finally we give a last alternative 
expression to g2 equivalent to (|12()|l at the present order: 

ff2(r) = (1 - 1/iV) p{Q)p{v) + {■.5p{v)6m 1)2 (121) 

where N is the mean total number of particles and p is the mean total density: 

p(r)-po(r) + (<5/5(r))3. (122) 



D. First order correlation function 



The first order correlation function of the field is defined as: 



31 (r) ^ (^t(r)^(o)) = (v/;5(0e'(^(°)-«(--»v^>- 



(123) 



As previously done we perform the calculation up to second order in the small parameters ei_2- We therefore expand 
\fp up to second order in 5p using H34(l . Note that we do not expand the exponential in 0{Q) — 9{r ), contrarily to 
what we did in the Hamiltonian: as r and are not neighboring points of the lattice anymore, the phase difference 
of the field can be arbitrarily large. The expansion in Sp gives rise to six terms: 



51 (r) 



Po 
1 



1/2, S 1/2,„N 

'^)Po (0) 



iAe\ 



-((5p(r)e 



iAdc- 



{5p\v)e 



iA0 



e'^^(5p"(0)-2,5p(r)e*^''(5p(0)) 



where we have introduced the following notations to simplify the writing: 

A0 = §{0) - 9{r) 
Sp{r) 



Po(r)' 



(124) 

(125) 
(126) 



We calculate the expectation values in this expression in two steps, first using the thermal equilibrium distribution 
for H2, and then including the corrections due to H^. 

The thermal expectation values corresponding to the quadratic Hamiltonian H2 are evaluated using Wick's theorem. 
One first expands the exponential in powers of A9, one calculates the expectation value of each term, and then one 
performs an exact resummation of the resulting series. This leads to the simple identities 



JAe\ 



{Sp'ir). 



-({Aef)2/2 



e-ii^'^rh/2(^Sp{r)tA9), 



iAe\ 



_ -{{Aef)2/2 



{6p^r)h + m{v)zA9)2y 



(Sp{r)e''^'Sp{0))2 = e-«^^)'>^/2 [(^Spir)Sp{0))2 + (<5/5(r)^A0)2(^A0^p(O))2] 



(127) 
(128) 
(129) 

(130) 



The expectation value of a product of the density fluctuation Sp and of the phase variation iA9 is particularly simple. 
In a classical field theory this expectation value would obviously vanish, as there is no crossed term in H2 between 
the density fluctuations and the phase. In the present quantum field theory this is not exactly the case as p and 9 
do not commute. To show that, we use the fact that the Bogoliubov mode functions Us(r), ^^(r) can be chosen to be 
real, so that dp and iA9 are linear combinations of 6^, 6j, P with real coefficients. As a consequence 



{Sp{r)iA9)2 = {5p{r)iA9)*2 = -{iA9Sp{r)). 



This leads to 



{Spir)^A9)2^-{mr),A9]h 



1 -5r 

2po{r)l 



D ■ 



(131) 



(132) 
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The same reasoning can be made for the other expectation value 

These expressions are second order in ei^2- An important consequence is that the product of such crossed phase-density 
expectation values in (|129I130|I are actually negligible at the present order of the calculation. The resulting form for 
gi, at the level of iJ2, is quite simple: 



l-i((A5p)2)2 + iscoria(r) 



(134) 



The notation A6p is similar to the one for the phase: 

ASj5 = 6p{0) - (5p(r) (135) 

and the 'scoria' comes from the crossed expectation value of A9 and 5p: 

1 1 
po(^ ^ MOW 



scoria(r) = (1 - ^r,o) ( —TXTd + „ /o^,n ) ■ (^^S) 



At this point we face the same apparent problem as in the calculation of g2. the 'scoria' scales as l/l^ and gives 
the wrong impression that our expression for gi will depend dramatically on / in the continuous limit 1^0. As in 
the case of 52, we solve this problem by expressing 6p and 9 in terms of the field A and putting the operators A, 
A''^ in normal order. We use (|118|l for the expression of Sp. For the difference of two phase operators, Q and the 
antihermitian operator a cancel so that 



where we have introduced the notations 



A6I = - f AA - AA^ 
2i V 



A(r 



(137) 



AW ^ ^ (138) 
Po (r) 

AA = A(0)-A(r). (139) 

After some calculations we arrive at 

{{ASpf)2 = (:(A5p)2:)2+scoria(r) (140) 

{{Ae)^)2 = (:(A0)2:)2 + iscoria(r). (141) 

As the 'scoria' is second order in ei_2, the exponential function of it can be expanded to first order. We then find as 
expected that all the l/l^ terms exactly cancel: 



l-li-.iASpf 



(142) 



The last step is to include the first order correction to gi coming from the cubic Hamiltonian H^. One then has 
to calculate expectation values with the thermal equilibrium density operator exp[—f3{H2 + H3)] to first order in H3. 
This thermal density operator can be viewed as the evolution operator during the imaginary time —ih(3 so that one 
can use first order time dependent perturbation theory to get 



(0)3 ~ (0)2 - / dr {e^"-H3e-^"-0)2 (143) 



where O is an arbitrary operator of the gas and where we have used the fact that H3 has a vanishing expectation 
value in the thermal equilibrium state for H2. One is back to the calculation of expectation values of some operators 
in the thermal state corresponding to H2. Wick's theorem can be applied. The resulting calculations are very similar 
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to the ones leading to (|142|l . but more involved and are detailed in the appendix^ The same phenomena takes place, 
that terms of a higher order than the present calculation can be neglected. One then gets 



iAe\ 



(144) 
(145) 



The first terms in the right-hand side of (|144|l and H145|l already appeared at the level of H2, the second terms are 
corrections due to H3 that we now take into account. There is no need to include H3 corrections to the other terms 
of p24|l since they are quadratic in 6p and are therefore already of second order. The expectation values of the phase 
9 and of the density fluctuations 6p have been calculated in subsection III Dl It was found that the expectation value 
of the phase operator is space independent so that {A.9)3 vanishes. The expectation value of the density fluctuations 
including the effect of H3 was given in (|95|l and is in general different from zero. Remarkably the whole effect on the 
correlation function gi of the first order correction in H3 is to replace po(j) by the total mean density p{r) defined in 

We write our final expression for the first order correlation function of the field, calculated consistently up to 2' 



gi{r) = y/p{r)p{0)exp 



(146) 



Note that we have inserted the contribution of the density fiuctuations inside the exponential factor, which is allowed 
at the order of the present calculation since this contribution is of order 2- 

What happens in the regime where a true condensate is present ? Both phase and density fluctuations are small, 
so that the exponential function in 1)1461) can be expanded to first order. We then express A9 and Ap in terms of the 
operator A and of the operator P. Since the Bogoliubov theory is usually considered in the canonical ensemble we 
remove the terms corresponding to the grand canonical fluctuations of the particle number. We then recover exactly 
the Bogoliubov prediction: 



9i 



Bog 



(r) = *,(r)vI/,(0) + (At(r)A(0)) 



(147) 



where 5'c(r) = y/Notpoir) + x(r)/-\/lVo is the condensate fleld. Amazingly the general formula 1)146(1 for quasi- 
condensates can be related to the Bogoliubov formula in the following very simple way: 



5i 



(r) = v/p(r)p(0) exp 



- 1 



(148) 



IV. EXPLICIT RESULTS FOR THE SPATIALLY HOMOGENEOUS CASE 

In this section we apply our approach to a spatially homogeneous Bose gas. The quasi-condensate density is then 
uniform: 



Po(r) 



90 



(149) 



The Bogoliubov equations (|50|) can then be exactly solved for any dimension of space and lead to Uk{r) = Uk e^^"^ jL^I"^ 
and Ufc(r) — Vk e'-^^ /L^/'^ with: 



Uk - Vk 



2m 



The corresponding eigenenergies are given by: 



1/4 



and Uk + Vk 



(^ + 2^^) 



2m 2m 



1/2 



2m 



2m. I 



1/4 



(150) 



(151) 
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A. Equation of state 



From the general expression (|103|l for the chemical potential of the gas we arrive in the thermodynamical limit at 



fJ' = P9o + 9Q 7:^-^[(wfe + Wfc) "fc + Wfe(wfc + Wfc)] (152) 
Jv (^■"■j 

where nj, — l/(exp(/3efe) — 1) is the mean occupation number of the Bogoliubov mode k. I? = [—tt/I, is the square 

domain of integration in the k space. The integral over the wavevector k does not contain any infrared divergence for 
any dimension of space. However the long wavevector behavior given by 

Vk(Uk + Vk) ^ (153) 
n 

gives rise to an integral convergent in ID and divergent in 2D and 31? in the 1^0 limit. This gives the impression 
that the result depends strongly on The solution of this paradox comes from the link between the bare coupling 
constant of the model potential in the discretized space and the low energy two-body scattering properties of the 
exact potential in the continuous space. This gives to go in two and three dimensions a dependence in / so that our 
expression for /i does not depend on I anymore in the I limit. In one dimension, the bare coupling go is simply 
equal to the actual coupling strength g for I — > and there is no divergence. A T = 0, the equation H152|l leads to: 

,^gp(l-^) (154) 

where ^ is the healing length defined in (|23|l . This agrees with the result of Lieb-Liniger in the weak interaction limit 
[28l |. In three dimensions, we refer to the appendix of '2^1 where the calculation has been done. One finds: 

3" = r-k — 

1-5 



g is the usual 3-D coupling strength given by : 



m 

where a is the exact potential scattering length. A more explicit form of H155|l is 



(156) 



where K = 2.442 .... It has to be noted that the difference between go and g is still small in the validity domain of 
our approach since it is a second order correction in ei^2- taking I ~ ^ one finds a/l ^ ^/ pl^- Replacing the first factor 
(7o in H152|) with the formula H155|) expanded up to second order in ei.2 gives: 

P = P9 + 9o[ j^yT ( iuk + Vkynk+Vk{uk + Vk) + ^jjt] (158) 
One can then safely take the I limit. At T = 0, the integration gives: 

^ = gp(l + ^^] (159) 



which is the same result as Lee and Yang's |33l | . In two dimensions, the low energy two-body scattering of a general 
short range potential is described by a single length a also named the scattering length. In a continuous space, the T 
matrix can be calculated in the low energy limit: 

{■k\T{E + ir])\k') ~ , , ^^""^^ ^ (160) 

^ ' ^ ' m[ln(2|a) +C-i7r/2] ^ ' 

where C — .57721 ... is the Euler constant, a is the scattering length, E — fi^k'^/m and r] 0+. We can also calculate 
the T matrix for the discrete delta potential defined by Eq. (|77|l which can also be expressed as: 

1/= |J|r = 0)(r = 0|. (161) 
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The general scattering theory gives the relations between the T matrix, the propagator G and the free propagator 
Go: 

T ^ V + VGV (162) 
G = Go + GqVG (163) 

Using these relations and Eq. 1)161(1 for the potential, we find: 

,m.AE + .,)|k') . i_,„;,.o|Gr(i; + i„)|r.O) - <""" 

The only term we need to calculate is the free propagator taken at the origin which is conveniently performed with a 
Fourier transform: 



(r = 0|Go(£; + z7?)|r = 0)= / „ ^ ^ . \,,^. 



(165) 



We split the square V into a disk of radius tt/1 and the complementary domain. Integration over the complementary 
domain gives simply a constant term in the low energy limit E <C h^/mP: 

27rh^ f dk 1 1 f dk 2G 



J= / TlTT^ 9 --7^ TT = ln(2) (166) 

where G — .91596 is the Catalan constant. The disk integration is straightforward and leads to the following expression 
for the T matrix: 

(k|T,.,(i. + z,)|k') ^ ^ ^ ■ (167) 

In — + — ^ J 



go 2TTn^ \ n J Ah^ 2TTn^ 



We now take Tgrid = T, where T is approximated by l(160() . in order to reproduce the low-energy scattering properties 
of the exact potential. This leads to: 

' - U±)-G^'-^). (168) 



50 27r?i^ V V Tray ' t: J 

Note that the condition (|26|) has to be satisfied in our approach. In two dimensions, this gives ti^ /mgo 3> 1 or using 
Eq. ((TC5|l : 

^lnf-)»l. (169) 
2n \a J 

We now show that the logarithmic dependence in / appearing in go, Eq. H168|) . exactly cancels the one appearing in 
the equation of state. Eq. H152|l can be rewritten as: 

^ ^ [iuk + Vkyrik + Vkiuk+Vk)] . (170) 



go Jv (27i")2 

In the thermal part, one can immediately take the I — > limit. In order to calculate the integral corresponding to 
the T = case, we use the same technique as for the calculation of go: the integration is made on a disk domain of 
radius n/l and we keep as a correction the integration over the complementary domain. The complementary domain 
integration is made by using the high wavevector behavior of Vk{uk + Vk), Eq. H153|) . This leads to: 



dk _ , _ _ , m/i 

;Vk{Uk + Vk) 









-1-2J 





Iv (27r)2 47rn 
Using Eas. (|168|) and (|171|l in Eq. (|170|) . we arrive at an implicit equation of state: 



(171) 



f 411^ \ f dk 



>- . t2 1'^ ^ — 2C+T ~ / yr^^^'^ + (i''^) 

47r?i \a^miie^'~'+^ ) J (2t:Y 

Remarkably this is identical to the result (20.45) obtained by the functional integral method in j7|i. A T = 0, one can 
show from the condition p^^ ^ 1, see Ea. (|26ll . that the validity condition of our approach is \n.{l/ pa?) 3> 47r. If one 
approximately inverts Ea. (|172|l neglecting constant terms and ln(ln(l/pa^)) with respect to \ti{\/ pa?), one recovers 
Schick's formula 1341. 
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B. Are density and gradient-of-phase fluctuations small ? 



As mentioned in section lll Bl our approach relies in particular on two assumptions: the assumption that the relative 
density fluctuation ei is small, and the assumption that the phase variation £2 between two neighboring points of the 
grid is small. 

Let us consider first the relative density fluctuations. Thanks to ()119|l their mean square value can be separated in 
two parts: 



''-^^^7^^ 7, ^^^^ 



where we have neglected with respect to l/l^ in the thermodynamical limit. The second term in p7;-{|l . involving 

the normal order, is expressed in terms of the Uk,Vk in the thermodynamical limit as 

(:(5/52(0):)2 2 f dk _ . . 

- ' [{uk+Vkfnk + Vk{uk+Vk)\ (174) 



pI Po Jv (27r) 

where the integration domain is I? = [—tt/I, tt/1]-^ . At zero temperature one introduces the change of variable q = 
in the integral: one finds that (|174|l is of the order of 1/po? in ID, of the order of ln(^/^)/poC^ in 2D and of the order 
of in 3D. Since I < ^ the second term in (|173|l is dominated by the first term and one has indeed 

(175) 



D • 



Pol 

At finite temperature we have to calculate the thermal contribution to Ij 174(1 involving the occupation number Uk ■ 

At a temperature kgT < fi we use the low momentum expansion of the Uk + Vk a-nd ej. and we find that the thermal 
contribution is {kBT/p,)-^^^{l/^)^ times smaller than 1/pol^. 

At a temperature ksT > that is A < ^, the treatment depends on the dimension of space. In ID the main 
contribution to the integral comes from the domain e^, ~ p, over which one can approximate the Bose formula by its 
low energy limit ksT/tk- This leads to a normal ordered fluctuation 1(174(1 of the order of ksT / {ppf)S_). This is larger 
than 1/ pqX so that the condition / < A then no longer implies that the first term 1/po^ in 1(173(1 is the dominant one. 
For convenience one can however adjust I to a value such that 



1 fcsT 1 



PqI M Po^ 
The condition for weak density fluctuations then becomes 

2 ksT 1 



(176) 



P PoC 



< 1. (177) 



Using pq ^ p and p gp we recover a condition already obtained in j5|| with a pure classical field approach. Note 
that this condition can be rewritten as ^ <C ?c where the coherence length of the field will be defined in ((187(1 . In 2D 
both the low energy domain efe < fc^T and the high energy domain > fc^T have an important contribution. In 
the low energy domain we approximate the Bose law by its low energy limit. In the high energy domain we keep the 
full Bose law but being then larger than /x, we approximate Uk + Vk by unity and by h'^k^/2m. This leads to a 
normal ordered fluctuation ((174(1 of the order of ln{kBT/ p)kBT/ {ppoS^'^), a quantity which is larger than 1/poA^. As 
in ID we therefore adjust I so that 

2 1 keT (kBT\ 1 

Ci In —-. (178) 

Pol M \ P J PoC 

In 3D the high energy domain > fc^T gives the dominant contribution so that the normal ordered expectation 
value 1(174(1 scales as 1/poA"^. This is dominated by the first term in ((173(1 so that the estimate ((175(1 applies as soon 
as Z < A, ^. 

Let us consider now the condition that the mean squared phase change over a grid cell: 

el ^ {{IVefh = TT- / 7^k\uk - Vkfink + 1/2) (179) 
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is much smaller than unity. The presence of the factor fc^ inside the integral, coming from the action of V, has the 
consequence that the contribution to the integral is dominated by the high energy domain. At zero temperature one 
can replace Uk — Vk by unity since the integral is dominated by wavevectors k ^ l/l > 1/^. This leads to 



el ^ -L. (180) 

whatever the dimension D. 

At a temperature fcsT < /i we estimate the thermal contribution by replacing Uk — Vk and by their low momentum 
approximations: the thermal contribution is then (Z/A)^+^(^/A)^ times smaller than the zero temperature result (|180|) 
and is therefore negligible since I < ^ < X. 

At a temperature fcsT > fi we use the high energy approximation replacing Uk — Vk by unity and Ck by h^k"^ /2m. 
Note that this works even in ID thanks to the presence of the fc^ factor in the integral (|179|l . This leads to a thermal 
contribution which is [l/XY'^^ times smaller than the zero temperature contribution l|180(l . and which is negligible 
since / < A. 

We conclude that the small parameter 62 of the theory, ensuring that there is a weak phase variation over a grid 
cell, is always given by (|180|l provided that the conditions (|22f) . H25|) and H26|l are satisfied. 

One may wonder if the condition that the corrections of the mean density due to the interaction between 
the Bogoliubov modes lead to an extra validity condition of our treatment. For the considered case of a spatially 
homogeneous gas it turns out that the answer to this question is no. One has indeed the remarkable identity in the 
thermodynamical limit: 

-{5p{v)h^-^{:5p{vf-h. (181) 

If the relative density fluctuations are weak, the relative correction to the density will be also weak. 

To end this subsection we discuss briefly the second order correlation function of the field (72 (r) • Restricting the 
general formula H121|) to the spatially homogeneous case in the thermodynamical limit we obtain 

52(r) = + 2p / — — [(ufe + Vkfuk + Vk{uk + Vk)] cos(k • r). (182) 

Limiting cases of this general formula can be compared to existing results in the literature. At zero temperature for 
a ID Bose gas one gets for r = 0: 

g^{Q)=p^(l-^\. (183) 



7rp^ 

This formula can be checked from p^ : the mean interaction energy per particle v is equal to 32(0) multiplied by 
g/2p, and v can be calculated in the weakly interacting regime by combining (3.29) of |2^ (relating v to the derivative 
of the ground state energy with respect to g) and (4.2) of |23 (giving the ground state energy in the Bogoliubov 
approximation). This exactly leads to H183|l . This prediction for 512(0) also appears in [35l |. 



C. First order correlation function 



Thanks to the general formula H148(l the first order correlation function of the field for the quasi-condensate is 
immediately related to the one of the Bogoliubov theory, here in the thermodynamical limit: 

ln[5i(r)/p] ^ - 1 = -i I [i^ul + vl) Uk + vl] (1 - cosk • r). (184) 

We have also taken here the continuous limit I ~* 0, which does not lead to any divergence. 

We concentrate our analysis to the ID case and we make the link with existing results in the literature. These 
existing results deal with the asymptotic behavior of gi for large r, where r is the absolute value of the spatial 
coordinate. At zero temperature, we find for r ^ ^: 

gi{r)^p[^) (185) 
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with ri = e ^/4 ~ 1.037^ where C = 0.57721 ... is Euler's constant [3g. This reproduces a result obtained in a 
non explicit way by |37j . At a finite temperature, gi/ p is the exponential of an integral of the form 

'•+°° A(k) 



fc2 



-[1 — cos(fcr)] 



(186) 



where the function A{k) is a regular and even function of fc, therefore behaving quadratically with k around fc = 
[38l| . Writing A(fc) as [A{k) — A(0)] + A(0) and splitting accordingly the integral one obtains for r much larger than 
both f and A: 



ln[5i(r)/p] = - + if + o(l/r") 
where the coherence length l^. ~ /cA^/tt coincides with the one of [s^l and the constant K is given by 



K = 



+ 00 



A{k) - A(0) 
P ■ 



(187) 



(188) 



Since A{k) is even one can show by repeated i nteg ration by parts that the remainder in H187|) tends to faster than 
any power law, contrarily to what is stated in [39|. 

Of course our formula gives access to gi for any value of the distance. This is illustrated in figure ^ where we 
have plotted the logarithm of gi as function of r/^ for various temperatures. As a consequence we can for example 




FIG. 1: First order correlation function of the field gi[z) for a repulsive ID Bose gas in the thermodynamical limit. The different 
curves correspond to various ratios of the temperature to the chemical potential: ksT / ^ = (solid lines), 1/15 (dot-dashed 
lines), 1/8 (dashed lines), 1/4 (dotted lines). We plot the logarithm of gi{z) multiplied by the parameter , where p is the ID 
spatial density and ^ = fi/y/rnfi is the healing length, so that we obtain a quantity depending only on z/£^ and ksT/fi in the 
weakly interacting limit. 



calculate the momentum distribution of the atoms: 



r+oo 

n(p) = 2 / dr gi{r) cos{pr /h) 
Jo 



(189) 



normalized here as / dpll{p) = 2'Khp so that n(p) is dimensionless. This is illustrated in figure HI where we have 
plotted the momentum distribution for various temperatures and for = 10. Using integration by part we can show 
that the behavior of 11 for large p is related to the fact that the third order derivative of (71 in r = 0"*" does not vanish: 



n(p)~?^!M^i^ with gf\Q+)^p''m^/{2h^). 



(190) 
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This prediction, valid at zero or finite temperature, agrees with the weak interaction limit of a recently obtained exact 
result based on the Bethe ansatz pX]| . At zero temperature we find that the momentum distribution diverges in p = 
as 

U{p)^^{np/hr (191) 

where ly = 1/{2ttp^) < 1. 




0.05 0.1 0.15 0.2 

FIG. 2: Momentum distribution of a repulsive ID Bose gas in the thermodynamical limit. n(p) is normalized as f dpn{p) = 
2Trhp where p is the ID spatial density so that n(p) is dimensionless. We plot n(p) as a function of p£,/h for various ratios of 
the temperature to the chemical potential; ksT / p, = 1/3 (dot-dashed lines), 1 (dashed lines), 10/7 (dotted lines). We have 
taken = 10 3> 1 where ^ — %/ ^mp is the healing length. The solid line is the large p limit: [Ti/p^)'^. 



V. CONCLUSION 

We have studied the thermal equilibrium of weakly interacting degenerate Bose gases in the regime of weak density 
fluctuations, the so-called quasi-condensate regime. The method can be considered as a Bogoliubov method in the 
density-phase representation of the field operator. 

In a first step one discretizes the real space in cells of size Z: / is small enough that the macroscopic properties of 
the gas are not affected by the discretisation, and I is large enough that each cell contains on the average a large 
number of particles. The macroscopic occupation of each cell allows to give a precise definition of the phase operator, 
as pioneered by Girardeau 

In a second step one performs a systematic expansion of the full Hamiltonian in terms of two small parameters, 
the relative density fiuctuations inside a cell and the phase change over a grid cell. This procedure leads to an exact 
expansion of the observables of the gas in the regime of weak interactions and low density fiuctuations, in ID, 2D 
and 3D. In particular it is free of any ultraviolet or infrared divergences and exactly matches the usual Bogoliubov 
predictions when the gas contains a true Bose-Einstein condensate. 

As a first application of the general formalism, we have given in this paper formulas for the equation of state of the 
gas, the ground state energy, the first order and second order correlation functions of the field. We have applied these 
formulas to the spatially homogeneous case in ID, 2D and 3D, recovering in this way known results, but obtaining 
also to our knowledge new results, like the full position dependence of the first order correlation function of the field. 
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APPENDIX A: EXPANSION OF THE HAMILTONIAN 



As explained in III Bl we expand the Haniiltonian (|28|l up to third order in powers of the small parameters ei and 
£2 defined in This will produce terms iy("i'"2) of order t^^e^^ with ni + 712 < 3. The expansion of the 

potential energy part -ffpot defined in H3U|) is very simple as it only involves the operator giving the density. The only 
point is to realize that the term l/l^ is €\ times smaller than the zeroth order density po- This leads to 



H-(0,0) 
^pot 




[/(r)-/i- 


, 50 
1- yPo 


(Al) 


rr(l,0) 
-"pot 


= Y,l''5p[u{v)~^i^ 


- ffoPo] 


(A2) 


o-(2,0) 
-"pot 


^ 2 

r 


" r -2 Po " 




(A3) 


rr(3,0) 
-"pot 






(A4) 



r 



The expansion of the kinetic energy part H31|l is more complicated as it involves also the phase operator 9, which 
furthermore does not commute with 5p. An expression slightly more convenient than 1)31(1 can be given for the kinetic 
energy. Thanks to the periodic boundary conditions one can freely shift the summation variable in the term of 1(31(1 
involving p_ j , so that 



kin 



2ml'' 



2Vp 



(A5) 



The calculation to zeroth order in £2 can be done first easily: using the expansion 1(37(1 to zeroth order, we get from 
((31(1 to all orders in £1: 



„(<+oo.O) 
-'^kin 



2m 



(A6) 



It involves a function of p only that it is easily expanded in powers of £1 using ((34(1 . A simplification occurs after 
summation over the lattice, as the matrix A is symmetric for the considered periodic boundary conditions: 



uAv = ^^(Au) V 

r r 

where u and v are arbitrary functions on the lattice. This leads to 



rr(0,0) 
-"kin ^ 


2m 


E 

r 


rr(l,0) 
^kin - 


2m 


E 

r 


77-(2,0) _ 
-"kin ^ 


2m 


E 

r 


r7-(3,0) 
^kin - 


2m 


E 



5^ 

Po 



Sp" 



3/2 



4\/A) Vpo 4/9o 

I5p^ 1 5p^ Sp 

8;pAV^---37,A-^ 

The second order term of vanishing order in £1 is also immediately obtained: 



(0,2) 
kin 



2ml 



'oPo,+i 



(A7) 

(A8) 
(A9) 
(AlO) 
(All) 

(A12) 
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The last second order quantity to calculate is Hj^^^'^ , which is first order in ei and first order in 62. There are four terms, 
two involving and two being their hermitian conjugate. One can then collect the terms to form commutators: 



(1.1) 

kin 



1/2 



[0, 



E 



1/2 



Po 

P0.+] 



. Sp+j] 

1/21 



l_ I Po-+j 
2 V Po 



1/2 



[0 



+3 



(A13) 
(A14) 



where we have used the commutation relation of p and 9, see (|20|) . 

We collect all the second order c-number contributions to the Hamiltonian H m & single energy functional of the 
density profile of the quasi-condensate, 



4mF 



PO 



1/2 



Po 

PO.+J 



1/2- 



(A15) 



The technique used to calculate H^^^'^ can be extended to the calculation of H^^^\ There are now three terms and 
their hermitian conjugates. Two of these terms, when combined with their hermitian conjugate, form a commutator 
that is calculated according to 12()() . The third term and its hermitian conjugate involve the expression 



5pid+ - o)Sp+ - 5p+{e+ - e)6p = 6p[e+,5p+] ~ {^p, d\8p+ 

which is a sum of two commutators, easy to evaluate. This leads to 



rr{2,l) _ fl (5/5 / 1/2 1/2 1/2. -l/2\ 

^kin -8^1^-[Po -p, Ap, j. 



(A16) 



(A17) 



(1 2) 

To calculate H^^^ ' wc first evaluate 



o-(< + oo,2) ^ JT- 



(A18) 



and we expand to first order in Sp, which leads to a sum of terms that are not individually hermitian. We then use 
the commutation relation H20|) to produce hermitian terms, e.g. 



6pid+3 - ef = [9+, - 9)5p{9+, -9)- j^{9+, - 



(A19) 



The last term of the right hand side of this expression is antihermitian and does not contribute to the final result 

t2 /.1/2 



(1,2) 
kin 



r,i \ Po Po,+j / 



(A20) 



Finally U^^^^ and Hyui^ vanish as the odd order expansion of exp[i(^^+j — 0)] is antihermitian. 



r(0,3) 



APPENDIX B: CORRECTIONS TO THE EQUATIONS OF MOTION DUE TO H-s 



The Hamiltonian i?3 gives rise to quadratic corrections to the equations of motion for bp and 9. In this appendix, 
these corrections are calculated explicitly and the thermal average is taken over the equations of motion with the 
Hamiltonian H2+H3 for the linear part and the Hamiltonian H2 for the quadratic corrections. This allows to calculate 
the first correction to the mean density due to H3. 

The corrections to the equation of motion for the density fiuctuations are given by: 



hdtSp\Hs 



■E 



Po,+j 
Po 



1/2 



5p 



Po 
Po,+j 



1/2 



5p+j 



(Bl) 
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where {A, B} stands for the anticommutator AB + BA of two operators. When we take the average with the 
Hamiltonian H2, we use the exphcit modal expansion of Sp and 9 given by Ea. (|55|l . The operator Q disappears since 
Eq. (|B1|I involves only differences of 9. Terms with P also disappear since {P)2 = 0. The expectation value of the 
product 9{r)dp{r'), where is written without the Q and 6p is written without the P, is actually purely imaginary: as 
Us and Vs can be chosen to be real, 9* — —9s, see Eg . (15711 . Since dt{Sp) is real, all imaginary contributions to it have 
to cancel so that the correction to the motion of ((5/5) due to H3 finally vanish when we take the thermal average: 

hdt{6p)\H,=0 (B2) 

The corrections to the equation of motion for 9 are more involved: 

^9Ah. = ^ (-:^SpA ( + i^.p^A (V^) - |1a f ^ 



Fortunately we can use the linear equations of motion H46I47|I to significantly simplify the above equation of motion. 
We rewrite the first term in parentheses of Ea. l|B3|l as: 



/ (5/5 \ 5p 



-SpA 



{u-fi + 350P0) ( ^ ) + '^Vp^ridte 



The second term in parentheses of Ea. lB3|l gives, as y/po solves the Gross-Pitaevskii equation: 



(B4) 



Sp^A{y/^) = —^{U - p + gapo). (B5) 



The sum of Eas. (|B4| ). ljB5|) and the third term in parentheses of Ea. (|B3|) lead to: 



?i A \ I Sp^\ Sp^ ^^r^n.Sp 

I — + U~p + gopo)[^]-9o^-m0)^- B6 



To rewrite the fourth term in parentheses of Ea. ljB3p . it is convenient to use the following identity: 

Vp^(0+, - 9f + V7^(^-. - of + 29 [^p^{9+, - 0) + Va^(^-, - 9)) = ^/p^ASl, - 0^) + VP^i^l, - 

- P(^A(^^9^)-9^Ai^)) 

leading to: 

^P^{9+j - 9f + - ef - f {9^A^ ~ 29 A (Vp^^) + A [■/^9^)) . (B8) 

Using this equality, the Gross-Pitaevskii equation H41|l and the equation of motion (|47() , the fourth term in parentheses 
of Ea. ljB3|l can be written as: 

E - + VP^^e-, - of) = [-^A + U-p + .gopo) (V^^^) - M^^. (B9) 



The sixth (and last) term in parentheses of Ea. ljB3p can also be transformed using the Gross-Pitaevskii equation 141|l : 



(,„-'^Apr-.rA..-"^) ^- -^A + + 5^ . ,B10, 



This leads finally to a rewriting of the thermal average of Eq. l|B3ll as: 



> \ A(5/52)2 ^ 1 ^ (<5/3^)2 - PoA^ ndt{hph 

2^n{dt69)\H3 = -7^ A + U - p + gopo — ^ + VPo(^ )2 - ^.^ _ - go 



(Bll) 
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The last term of this expression can be calculated using Ea. (|56|l . The harmonic modes do not contribute since the 
expectation value of products of bs and is time independent. We are left with: 

dt{QP)2 = dt{Ql0)P + t^P')2 = ^(P')2 (B12) 



ndt{05p)2 „ ,/p2\ ^ ^ /T31q^ 



which gives: 



Po 

As a conclusion, the quadratic correction to the first equation of motion can be written as in (|60|l if one uses the 
identities: 



^i^B + B^y = ^(2B^B + B^+B^^ + ^y (B14) 



APPENDIX C: THE MEAN VALUE OF dtO VANISHES AT EQUILIBRIUM 

As the field degree of freedom Q, that is the global phase of the field, is not subject to a restoring force in H2, it 
is not totally obvious that the perturbation H3 cannot set it into permanent motion. We therefore check this point 
explicitly here. 

The first step is to calculate the mean value of P to first order in H^. We approximate the unnormalized density 
operator of the gas at thermal equilibrium to first order in H3 using perturbation theory: 

^ = g-/3(if.+ff3) = _ e-(^--)«^H3e-^^^ + . . . . (Cf ) 



P commutes with H2 and has a vanishing mean value in the thermal state corresponding to H2 so that, to first order 
in H3: 

(P)3 = -{I3PH^)2. (C2) 
The Hamiltonian H3, is a polynomial of degree 3 in P: 

H3 ^ Ao + AiP + A2P^ + A^P^ (C3) 
where the Ai are still operators with respect to the harmonic oscillator variables bg. This leads to 



{P)3 = ~(3 



(Ai)2(P^)2 + (^3>2(P">2 . (C4) 



From Wick's theorem, (P^)2 = 3(p2\2 



2- 



In a second step we calculate (dQ/dt) to first order in iJ; 



3- 



{dQ/dt) ~ {dp{H2+H3))3 

~ Mo(^>3 + (^l>2+3(A3>2(P')2 (C5) 

where the terms coming from dpH^ are calculated in the thermal state for H2 since they are already first order in 
the perturbation. From the value of (P)3 obtained from l|C4|l and from we obtain the desired result: 

(dQ/dt) 3 = (C6) 



to first order in H^. 
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APPENDIX D: AN EQUATION FOR {a, A} 

In this appendix, we derive the partial differential equation (|91() . We first note that _B„ , being a sum of eigenmodes 
of the operator CqPi obeys the differential equation for the evolution governed by H2'. 

^^dJfi]=Cop(fA (Dl) 



We project this equation orthogonally to 0o f^nd along (jjQ, so that we get the quantum analog of equations (E9), 
(ElO) of [2^, with the simplification that (/)o(r) is real: 



= 50/^000 (S„ + St ) = E (A + At) . (D3) 

r r 

We have introduced the projection matrix 

{r\Q\r')=5r,r'~l°Mr)Mr')- (D4) 



As a is antihermitian, the source term in (|D2p vanishes and one can replace i?„ by A in (|D3|I . 

We use these two equations of motion to calculate the first order time derivative of A{r) = {{a, A(r)})2. We do not 
give the intermediate result. As A is real here, we have the property 

({d,At(r)}>2 = -({a,A(r)})2. (D5) 

As A is orthogonal to (j)Q one has 

QA = A. (D6) 

All this leads to lEJ. 



APPENDIX E: INTERPRETATION OF x IN THE NUMBER CONSERVING BOGOLIUBOV 

APPROACH 

We assume here that the gas is a quasi-pure condensate so that 0o is now the condensate wavefunction in the Gross- 
Pitaevskii approximation. We then show that x(r)/A"o, where x is defined in (|94|l . essentially coincides with the lowest 
order deviation of the exact condensate wavefunction from the Gross-Pitaevskii prediction (pQ. This deviation was 
calculated in p^ . 

We split X in a component orthogonal to 00 and a component colinear to c/fo : 

X(r) =7(/)o(r)+x±(r). (El) 

The component 7 has a simple physical interpretation: we sum H95|l over r after multiplication by . If wc omit the 
grand canonical term (absent in the canonical treatment of 'l5]) we obtain 

N ^ Na + 2-/ + SN (E2) 

where 

<5iV = ?^E(A^A)2 (E3) 

r 

exactly coincides with the mean number of non-condensed particles predicted in . The physical interpretation of 
27 is then simple: 



6N0 = 27 



(E4) 
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is the correction to apply to the pure condensate prediction for the number of condensate particles in order to recover 
the correct Bogoliubov prediction. Applying to H94I) the matrix Q (|D4() projecting orthogonally to (j>o we obtain 



^2 

-—A + U + gopo - fi 
2m 



2Q9oPoX± + Q ( 2goPo70o + ^5 ) = 0. (E5) 



We modify slightly the form of the source term S, eliminating the anticommutator: 



{At(r'), A(r)} = 2At(r')A(r) + ^{r\Q\r'). (E6) 



This leads to the system 

with the effective source term 

= 9oPo{r)Mr){5No~l)+goNoMr) (2(At(r)A(r))2 + (A2(r))2)-Z^ ^.gopo(r')<^o(r')((A(r') + At(r')) A(r))2 

where we used the fact that here (A^)2 = (A^^)2 since the condensate wavefunction is real. Equation l|E7|l is the 
steady version of Eq. (95) of 0| , which gives N times the correction to the condensate wavefunction, and the source 
term (|E8|I exactly coincides with the one of Eq.(96) of if one realizes that N = Nq, so that SNq = —SN, in the 
systematic expansion used in . 

APPENDIX F: CORRECTIONS TO gi DUE TO THE CUBIC HAMILTONIAN 

We calculate the corrections to the first order correlation function due to using the perturbative formula (|143|l . 
A first remark is that 

i/3(T) = e^^-H^e-^'^' (Fl) 

is still cubic in the operators bg since one has 

e^^^ke"^"^- = e-^'^bs (F2) 

e^^'ble--'"' = e'-'^bl (F3) 

where is the energy of the Bogoliubov mode s. The second step is to use Wick's theorem to calculate the expectation 
values in the thermal state corresponding to Hamiltonian H2. One can derive the general formulas 



{AiA2A3e'^% = [{AiA2A3iA9)2 + {AitA9)2{A2iA9)2{A3tAe)2]e-^(^^'>'> '>-/^ (F4) 
(AiA2A3A4e*^«)2 = f{A,A2A3A4 



u 2 



l + l{{A9f) 



^{AiA2A3A4A9y) 



+ {A,iA9)2{A2tA9)2{A3tA9)2{A4tA9)2je-''(^''>'>-/^ (F5) 

where the Ai are linear in Sp and 9 and have a vanishing mean value. A last point is to realize that some of the 
obtained terms contain a larger number of factors equal to than other ones. Since A9 scales as 1/ \/poF^, see e.g. 
the expression of 9 in terms of the mode functions Us, Vs in (|57|l . the terms with an excess of A9 factors are higher 
order in the expansion and are therefore negligible. Note that strictly speaking this argument is correct provided that 
each factor {AiA9)2 remains bounded whatever the distance from to r. This can be checked indeed to be the case: 
from the form of H3 one sees that Ai is either dp or the phase difference between two neighboring points of the grid. 
One can therefore use the approximate identities: 

{AiA2A3e'^')2 ~ (AiA2A3iA0)2e-<(^«)'>^/2 (F6) 
{AiA2A3Aie'^')2 - (AiA2A3A4)2e-«^«)'>^/^. (F7) 
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This immediately leads to the identities H144I145|I . 
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